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Abstract 

A  distributed  parameter  model  describing  spatially-dependent  hepatic  processing  of  the 
chemical  compound  2,3,7,8-tetrachlorodibenzo-p-dioxin  (TCDD  or  dioxin)  has  previously 
been  reported  [1].  The  mathematical  system  consists  of  coupled  nonlinear  partial  and  ordi¬ 
nary  differential  equations  with  delays.  In  this  paper  we  investigate  the  qualitative  behavior 
of  the  system  over  a  six-hour  time  period  following  a  subcutaneous  injection.  A  brief  sum¬ 
mary  of  the  model  is  also  given. 

Keywords  and  phrases:  Distributed  parameter  models,  pharmacokinetic  models,  2, 3,7,8- 
tetrachlorodibenzo-p-dioxin,  liver  transport  models,  delay  equations. 


1  Introduction 

Numerous  physiologically-based  pharmacokinetic  (PBPK)  models  have  been  developed  to  de¬ 
scribe  the  uptake  and  elimination  of  the  environmental  toxin  2,3,7,8-tetrachlorodibenzo-p-dioxin 
(for  example,  see  [2] — [6] ) .  While  most  PBPK  models  assume  well-mixed  compartments,  recent 
work  [1,  6,  7]  has  focused  on  the  development  of  advanced  techniques  to  account  for  spatially- 
dependent  tissue  responses. 

In  rodents,  TCDD  has  been  shown  to  bind  to  two  hepatic  proteins:  the  aryl  hydrocarbon 
(Ah)  receptor  [8]  and  an  inducible  protein,  cytochrome  P4501A2  (CYP1A2)  [9,  10].  Following 
administration  of  TCDD,  binding  to  the  Ah  receptor  results  in  a  dose-dependent  induction  of 
CYP1A2  [9,  10];  therefore,  this  protein  is  present  at  both  a  basal  (non-induced)  and  induced 
level  in  the  presence  of  TCDD.  The  liver  response  to  dioxin  is  not  homogeneous,  resulting  in 
preferential  localization  around  centrilobular  regions  at  low  doses  [6].  Lumped  parameter  PBPK 
models,  which  assume  a  uniform  concentration  of  the  toxin  throughout  the  liver,  cannot  account 
for  this  known  heterogeneity  in  cellular  response. 

A  partial  differential  equation  model  to  describe  the  spatially-dependent  hepatic  processing 
of  TCDD  has  been  proposed  [1,  11]  and  is  briefly  summarized  in  Section  2.  In  Section  3  the 
numerical  scheme  used  in  solving  the  resultant  system  of  equations  is  given.  Our  preliminary 
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computational  results  on  the  qualitative  behavior  of  the  system  are  discussed  in  Section  4.  In 
particular,  we  examined  the  behavior  of  the  system  on  the  time  interval  from  zero  to  six  hours, 
prior  to  TCDD-induced  CYP1A2  synthesis.  We  also  studied  the  dependence  of  the  solution 
on  two  important  model  parameters,  the  dispersion  number,  T>n,  and  a  circulatory  delay,  rc. 
Finally,  in  Section  5  we  discuss  future  research  directions. 


2  The  TCDD  Model 


A  mathematical  model  (1)  has  been  developed  [1,  11]  to  describe  pharmacokinetic  and  pharma¬ 
codynamic  properties  of  TCDD.  A  convection-dispersion  equation  (la),  based  on  the  work  of 
Roberts  and  Rowland  [12],  characterizes  the  transport  of  blood  elements  in  the  liver  sinusoidal 
(blood)  region.  Throughout  this  discussion,  the  dimensionless  spatial  variable  x  takes  on  values 
in  the  range  [0,1];  x  =  0  corresponds  to  the  liver  inlet,  while  x  =  1  corresponds  to  the  out¬ 
let.  Uptake  of  dioxin  into  the  hepatic  cells,  called  hepatocytes,  is  assumed  to  occur  by  passive 
diffusion.  The  model  includes  the  dynamics  of  TCDD-binding  with  two  intracellular  hepatic  pro¬ 
teins,  the  Ah  receptor  (lc)-(ld)  and  an  inducible  microsomal  protein,  CYP1A2  (le)-(lf).  The 
induction  mechanism  is  described  in  terms  of  the  fractional  occupancy  of  the  Ah  receptor  at  a 
previous  time,  t  —  Tr,  to  account  for  the  many  intracellular  processes  which  must  occur  before 
an  increase  in  CYP1A2  concentration  is  realized.  Elimination  in  the  liver  (by  metabolism  and 
biliary  clearance)  is  assumed  to  be  a  first  order  process. 

A  well-mixed,  combined  venous/arterial  blood  compartment  ( lg) ,  which  includes  a  loss  due  to 
the  uptake  and  elimination  of  TCDD  in  the  rest  of  the  body,  completes  the  system.  A  circulatory 
lag,  rc,  accounts  for  the  time  delay  in  transport  of  blood  elements  from  the  exit  of  the  liver  to 
the  venous  measurement  location. 

The  mathematical  system  under  consideration  is  as  follows: 


( VB  +  VD  ^S.) 
JuD 

dCUH 

dt 
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CB(t,0)  =  Ca(t), 

QCB(t,  1)  -  QVN^f(t,  1)  =  Aq2(t), 
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where  C  =  [(  '/; •  CK ,, .  CA /.  /  .  C,\ /, ,  <Jpr  /.,  Cpr ,  C'„  7  .  In  Eqn.  (lh),  $  and  Aq2  are  assumed 
known. 

We  note  that  gAj l  and  gpr  are  saturating  nonlinearities  modified  from  the  usual  product 
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terms, 


(.</:  z)  =  k+1yz , 

gpr(y,z)  =  k+2yz ,  for  y,  z£t. 


(2) 

(3) 


arising  from  the  law  of  mass  action  in  chemical  kinetics.  Specifically,  we  assume  that  within  a 
certain  range  of  concentrations  the  system  behaves  according  to  the  nonlinearities  prescribed  by 
(2)  and  (3)  but  eventually  saturates;  i.e. ,  due  to  the  availability  of  binding  species,  we  assume 
the  rates  of  formation  of  Ah-TCDD  complex  and  CYP1A2-TCDD  complex  are  bounded. 

The  mathematical  model  (1)  consists  of  a  nonlinear  system  of  partial  differential  equations 
with  delays  and  questions  of  well-posedness  are  of  interest.  Detailed  general  theories  of  existence, 
uniqueness,  and  continuous  dependence  for  certain  nonlinear  parabolic  systems  (including  the 
TCDD  model)  are  given  in  [11,  13].  Furthermore,  a  general  result  for  related  least-squares 
parameter  estimation  problems  is  given  in  [14].  These  results  were  obtained  using  a  weak  or 
variational  formulation  of  the  problem  and  are  based  on  ideas  from  the  works  of  Banks  et  al.  [15]- 
[16]  for  nonlinear  hyperbolic  systems. 

The  summary  given  above,  while  brief,  is  included  to  provide  the  reader  with  a  general 
understanding  of  the  complex  dynamics  of  the  system  under  investigation.  The  reader  is  referred 
to  the  aforementioned  works  [1,  11,  13,  14]  for  complete  discussions. 

3  Numerical  Methods 

In  this  section  we  present  an  overview  of  the  numerical  methods  used  in  our  simulations.  All 
coefficient  matrices,  functions,  vectors,  etc.  are  as  described  in  [11], 

Our  numerical  scheme  is  based  on  a  weak  formulation  [13]  of  the  problem  given  in  Eqn.  (1). 
To  obtain  the  weak  form,  we  multiplied  the  ith  equation  by  a  function  ipi  in  a  ’’suitable”  class  of 
test  functions  and  then  integrated  in  space  in  the  first  six  equations,  followed  by  integration  by 
parts  in  Eqn.  (la)  only. 

For  ease  of  notation,  we  define 

y  =  [Cb,  Cuh,  CAh-T,  CAh,  Cpr-T ,  Cptt  Ca]T  =  [i/1,  ...  ,  t/7]T- 


3.1  Finite  Element  Formulation 

Let  0  =  xo  <  x\  <  . . .  <  xpf  =  1  be  a  uniform  partition  of  the  interval  [0,  1]  into  N  subintervals 
of  length  h  =  1/N.  We  take  as  basis  elements  the  piecewise  linear  continuous  functions,  <j>j, 
j  =  0, . .  . ,  N,  defined  by 

Xj-l  <  X  <  Xj, 

Xj  <  X  <  Xjp  l, 

0  <  x  <  Xj-i  or  Xj  +  i  <  x  <  1. 

We  define  the  Galerkin  finite  element  approximations  by 

.</:v  (':*•)  =  1 

for  i  —  2, . .  .,  6,  where  the  basis  elements,  <f>jt  are  as  described  above  and  yf  m  y{. 

Next  we  substitute  the  finite  element  approximations  (4)  into  the  weak  form  of  the  equations. 
Let  yN  (t)  6  JJ  JV+5(JV+1)+1  such  that 

yN(t)  =  [a1(t),a2(t),...,a6(t),y7{t)]T. 

The  finite  dimensional  system  we  obtain,  in  terms  of  the  time-dependent  coefficients  of  the 
Galerkin  approximations,  is  given  by 

MyN(t)  =  AyN (t)  +  GP{yN (i))  +  T{t)  +  ADyN (t  -  rc)  +Qs{yN{t  -  r,-)),  (5) 


<f>j  (*)  = 


h  5 

Vj. 1-1-37 
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where  the  matrices  M  and  A  are  elements  of||l,v+5(JV+1)+1)>'(',v+5(',v+1)+1)  and  the  vector-valued 
functions  Qp,  T,  Ad,  and  Qs  are  elements  ofj|jiv+5(JV+1)+1. 

The  initial  condition,  ,  for  the  semi-discrete  problem  (5)  is  taken  as  the  L 2  projection  of 
the  original  initial  condition,  <F,  onto  the  finite  element  space. 

3.2  General  Algorithm 

In  these  preliminary  investigations,  we  seek  solutions  on  the  time  interval  from  zero  to  six  hours. 
For  simplification,  we  have  assumed  that  no  TCDD  is  present  in  the  system  on  the  interval  [-6,0] 
hours.  By  making  this  assumption,  we  can  thereby  ignore  the  induction  delay  term,  Qs ,  for  the 
purpose  of  this  discussion. 

The  semi-discrete  problem  we  consider  is  given  by 

MyN (t)  =  AyN (t) +gp(yN (t))  +  F(t)  +  ADyN {t  ~  tc), 

yN{s)  =  vo(s),  se[-rc,o]. 

The  entries  in  the  coefficient  matrices  M  and  A,  consisting  of  certain  combinations  of  inner 
products  of  basis  elements  and  their  derivatives,  were  calculated  analytically  rather  than  through 
use  of  a  numerical  integration  scheme.  The  nonlinear  vector  function  Qp  was  treated  in  a  similar 
manner. 

The  “method  of  steps”  [17]  for  ordinary  differential  equations  was  used  to  solve  the  problem 
on  each  delay  interval  of  length  rc  (rc  «Tr).  The  general  algorithm  is  described  below: 

Algorithm. 

•  Set  T  =  final  time  (T  <  7>) 

•  Form  the  coefficient  matrices  M  and  A 

•  Solve  the  linear  system  defining  the  initial  condition  y in  the  finite  element  space 

•  Set  tO  =  0  and  tf  —  {  ^  T<[  " ^  ^ . 

[  1  ,  otherwise 

•  do  while  tO  <  T 

—  Solve  on  [70,  tf] 

M'i/N(t)  =  AyN (t)  +  gp(yN(t))  +  T{t)  +  AdVo  (*  -  Tc ) 
yNm  =  dm- 

-  Set  y$  (t)  =  yN  ( t ) 

-  to  =  to  +  Tc 

I J  /  1 0  +  Tc  ,  to  +  Tc  <  T 
J  T,  otherwise 

•  end 

In  order  to  evaluate  y^ ,  we  stored  the  computed  solution  throughout  the  previous  delay 
interval  and  used  an  interpolation  routine  to  find  the  value  of  the  solution  at  time  t  —  rc.  In  fact, 
only  (y^)jv  and  (y^  jgjv+e  were  stored,  as  these  are  the  only  data  required  for  the  linear  delay 
term,  Ad- 
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Figure  1:  Numerical  solution  obtained  for  Cs(f,  1)  varying  the  number  of  basis  elements,  Nx, 
with  T>n  =  0.5.  The  solutions  with  Nx  =  16  and  32  are  identical  on  this  scale,  indicating 
convergence  of  the  Galerkin  approximations  with  16  basis  elements. 

3.3  Implementation 

The  computer  code  was  written  in  MATLAB  version  5.1  (The  Mat  1 1 Works,  Inc.,  Natick,  MA) 
and  computations  were  carried  out  on  a  Sun  Sparc  Ultra  II  workstation  and  a  personal  computer 
with  a  180  MHz  Intel  Pentium  Pro  processor.  The  MATLAB  routine  ’odel5s’  [18]  was  used  for 
time  stepping,  which  is  a  variable  order,  variable  step  method  based  on  numerical  differentiation 
formulas.  The  relative  and  absolute  error  tolerances  were  set  to  1  x  10-b.  Since  this  is  a  variable 
step  method,  at  time  t  the  solution  over  the  previous  delay  interval  had  to  be  interpolated  in 
order  to  determine  the  value  of  the  solution  at  time  t  —  rc.  MATLAB  s  interpolation  routine 
’interpl’  was  used  with  the  ’spline’  option.  This  method  fits  the  data  with  cubic  splines  between 
data  points,  so  that  each  segment  of  the  curve  fit  has  at  least  the  same  first  and  second  derivatives 
as  the  ones  adjoining  it.  The  maximum  order  of  the  integrator  was  set  to  two  in  order  for  it  to  be 
in  the  range  of  accuracy  of  both  the  interpolation  routine  and  the  finite  element  approximations. 

3.3.1  Convergence  of  the  Numerical  Scheme 

The  theoretical  results  presented  in  [11,  14]  guarantee  convergence  of  the  Galerkin  finite  element 
approximation  scheme.  However,  we  must  determine  the  number,  N ,  of  basis  elements  to  be 
used  in  meaningful  simulations.  We  computed  solutions  over  the  time  span  from  0  to  1  hour 
with  fixed  Vn  and  rc,  and  varied  N .  Numerical  results  obtained  using  a  small  number  of  basis 
elements  and  t  very  close  to  zero  exhibited  dynamics  which  were  only  eventually  resolved  by 
decreasing  the  size  of  the  spatial  grid.  This  behavior  is  shown  in  Figure  1  with  Vn  =  0.5.  For 
values  of  the  dispersion  number  of  interest  (0.5-100),  we  found  that  increasing  Zf/v  generally 
resulted  in  a  decrease  in  the  number  of  basis  elements  required.  For  example,  convergence  was 
obtained  with  16  basis  elements  for  the  case  V jy  =  5  (Figure  2),  as  compared  to  just  8  basis 
elements  for  V jv  =  10  (Figure  3). 

Increasing  the  dispersion  number  has  the  effect  of  smoothing  concentration  gradients  along 
the  cylinder’s  length  [12];  that  is,  spatial  variations  in  concentrations  decrease  with  increasing 
Vj\r,  and  therefore  a  less-refined  spatial  grid  is  required  to  capture  system  dynamics. 
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Figure  2:  Numerical  solution  obtained  for  Cs(t,  1)  varying  the  number  of  basis  elements,  Nx, 
with  T>n  =  5.  The  solutions  with  Nx  =  16  and  Nx  =  32  are  identical  on  this  scale,  indicating 
convergence  of  the  Galerkin  approximations  with  16  basis  elements. 


Figure  3:  Numerical  solution  for  Cs{t,  1)  varying  the  number  of  basis  elements,  Nx,  with  T>n  = 
10.  The  solutions  with  Nx  =  8,  16,  and  32  are  identical  on  this  scale,  indicating  convergence  of 
the  Galerkin  approximations  with  8  basis  elements. 
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Model  Simulations 


4.1  Parameter  Estimates,  Boundary  and  Initial 
Conditions 

The  majority  of  parameters  used  in  these  simulations  were  obtained  from  the  literature  (see 
[11]).  However,  in  order  to  qualitatively  estimate  the  time-dependent  exit  flux  term,  Aq2{t), 
we  implemented  a  numerical  scheme  for  the  TCDD  model  satisfying  a  homogeneous  Neumann 
boundary  condition,  1)  =  0.  From  a  balance  of  fluxes  we  have 

Q(CB(t,  1)  -  1))  =  W>(t):  (6) 

from  which  we  obtained  a  discrete  approximation  to  Aq2: 

Aq2{ti)  =  QCB(ti,  1). 

This  data  was  then  interpolated  and  used  as  an  estimate  for  Aq2  when  simulating  model  behavior 
with  the  Robin  boundary  condition  (6). 

All  simulations  were  carried  out  assuming  zero  initial  conditions  for  concentrations  involving 
TCDD.  Constant  values  were  assumed  for  the  initial  levels  of  the  Ah  receptor  and  the  induced 
binding  species,  CYP1A2.  Uptake  of  TCDD  into  the  arterial/ venous  blood  supply,  I(t),  was 
assumed  to  be  from  a  subcutaneous  dose  of  300  ng/kg  body  weight.  This  was  described  as  a 
first-order  process  with  rate  constant  0.05/hour,  as  in  the  work  of  Anderson  et  al.  [5], 

The  qualitative  behavior  of  the  system  was  investigated  while  values  for  the  axial  disper¬ 
sion  number,  Xhv,  and  the  circulatory  delay,  rc,  were  varied,  holding  all  other  biological  and 
physiological  parameters  fixed. 

4.2  Predicted  Behavior 

The  system  was  studied  for  Vn  —  50,  rc  =  1  minute,  and  0  <  t  <  6  hours.  The  most  interesting 
behavior  occurs  within  the  first  two  hours.  As  the  concentration  of  TCDD  in  the  liver  blood 
increases  (Figure  4),  free  dioxin  in  the  cells  is  quickly  bound  to  the  Ah  receptor  and  CYP1A2. 
The  Ah  receptor  quickly  becomes  saturated  (Figure  5)  due  to  its  low  binding  capacity.  The  con¬ 
centration  of  the  AhR-TCDD  complex  continues  to  rise,  however,  as  newly  synthesized  protein 
binds  any  available  TCDD.  Inactivation  of  the  Ah  receptor  is  modeled  as  a  first-order  process 
(concentration  dependent),  so  that  at  low  concentrations  the  zero-order  process  of  protein  syn¬ 
thesis  dominates;  that  is,  there  is  a  net  increase  in  cellular  levels  of  receptor  available  for  binding. 

The  second  dioxin-binding  species,  CYP1A2,  takes  longer  to  become  saturated  than  the 
Ah  receptor  because  of  its  higher  binding  capacity  (Figure  6).  After  roughly  two  hours,  the 
concentration  of  free  CYP1A2  is  essentially  depleted,  as  any  newly  synthesized  protein  rapidly 
binds  to  free  dioxin  in  the  cells.  A  comparison  of  unbound  TCDD  in  the  liver  blood  to  that  in 
the  hepatocytes  is  given  in  Figure  7.  The  binding  of  TCDD  to  intracellular  proteins  causes  an 
initial  decrease  in  the  cellular  concentration  of  free  dioxin  relative  to  that  in  the  blood.  Once 
the  two  proteins  become  saturated,  however,  the  unbound  concentration  of  TCDD  in  the  cells 
begins  to  rise,  but  at  no  time  does  the  cellular  concentration  exceed  that  in  the  blood.  In  other 
words,  for  this  set  of  model  parameters  the  hepatocytes  act  as  a  sink. 

4.3  Dependence  on  the  Circulatory  Delay 

The  rate  of  change  of  the  combined  arterial/venous  blood  compartment  is  given  by  the  differential 
equation 

-  TC,  1)  -  Ca(t))  +  I(t)  -  keCa(t),  (7) 

with  initial  condition 

0.(0)  =  0.  (8) 
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Free  CYP1 A2  versus  CYP1 A2-TCDD  Complex:  Nx  =  16,  xc  =  0.016667  hrs,  DN  =  50 


Figure  6:  Concentrations  of  free  CYPIA'2  and  the  CYP1A2-TCDD  complex.  The  protein  be¬ 
comes  saturated  after  approximately  two  hours. 


Free  TCDD  in  liver  blood  versus  hepatocytes:  Concentrations  for  N  =  16,  =  0.016667  hrs,  DN  =  50 


Figure  7:  Concentrations  of  free  TCDD  in  the  liver  blood  ( CUB )  and  hepatocytes  (CUH)  over  a 
six  hour  period. 
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If  we  assume  a  constant  rate  of  uptake  of  TCDD  from  a  subcutaneous  depot  over  a  finite  time 
interval  [0,7],,], 

Cin,  t  <  Tin 


m  = 


0,  t  >  Tin, 


Eqn.  (7)  can  be  rewritten 

dCa(t) 
dt 

where  ci  =  ^r-,  fi  =  c\  +  ke,  and 

X[a,b]^  =  {  o| 

We  consider  the  only  two  possibilities  of  interest,  rc  =  0  and  rc  >  0. 


—  Cl  Tc ,  1)  fi Ca{ /)  H-  (^) 5 


t  e  fab], 

otherwise. 


(9) 


(10) 


Ignoring  the  Circulatory  Delay 

For  the  case  rc  =  0,  the  analytic  solution  to  (9)  is 

r/  m  =  /  “  e“^  + Cl  to  1)*,  0  <  t  <  Tin 

From  (10),  it  is  clear  that  for  any  t  >  0,  the  concentration  of  TCDD  at  the  exit  from  the  liver, 
C'ij  (t .  1),  contributes  immediately  to  the  predicted  arterial  blood  concentration  Ca(t). 

4.3.1  Predicted  Solution  with  Circulatory  Delay 

We  now  consider  the  case  rc  >  0.  Let  t  £  [0,T],  with  T  <g  Tn.  We  assume  a  zero  background 
concentration  of  TCDD  in  the  system;  in  particular,  Cb(1,  1)  =  0  on  [—  rc,0].  On  the  initial 
interval  of  integration  [0,  rc],  the  rate  of  change  of  Ca  is  given  by 

^=jf-{t)  =  ~nCa{t)  +  Cin\[0,T,n)(t), 

Ca{  0)  =0. 


The  analytic  solution  is 


C„(l)  -  —(\  <  »').  I  fc  [0,r(.]. 

ft 


On  the  next  interval  of  integration  [rc,2rc],  the  governing  differential  equation  for  Ca  is  given 
by  (9)  with  the  initial  condition 


The  solution,  for  t  £  [rc,  2rc],  is 


Ca(rc )  =  — (1  -  e~ 
ft 


Ca(t)  =  e-^-^Ca(Tc)  +  ^(l-e-^-^) 

ft 

+ci  f  e“^‘“s)C'B(s  -  rc,  l)ds 

*  rc 

=  g^(t^)Ca(Tc)  +  Ca{t  -  tc)  +  Cl  [  e-^-^CB(s  -  rc,  1  )ds. 

The  contribution  to  Ca  from  the  liver,  CB{t,  1),  does  not  appear  until  t  >  rc  >  0.  As  shown 
in  Figure  8,  the  concentration-time  profile  of  the  combined  arterial/ venous  blood  compartment 
illustrates  this  phenomenon  and  its  impact  on  the  solution  over  a  period  of  ten  delay  intervals. 
The  solution  follows  an  apparent  exponential  curve  on  [0,rJ;  after  this  time,  a  marked  increase 
in  concentration  is  seen  as  the  initial  input  of  dioxin  to  the  liver  (that  which  has  not  been  taken 
up  by  the  hepatocytes)  finally  exits  the  organ  and  returns  to  the  general  circulation.  As  time 
progresses  this  sharp  demarcation  in  solutions  at  the  end  of  each  delay  interval  diminishes. 


/I  Tc 


)• 
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C  (t)  for  Dm  =  5,  X  =  0.01 6667  hrs,  N  =  1 6 

a'  '  N  c  x 


Figure  8:  Predicted  concentrations  of  TCDD  in  the  combined  arterial/venous  blood  supply,  Ca  , 
for  a  period  of  ten  circulatory  delay  intervals. 


4.3.2  Predictions  of  the  Model  with  and  without  Circulatory  Delay 

The  concentration  of  TCDD  in  the  arterial/venous  blood  supply  ( Ca )  is  illustrated  in  Figure  9 
for  different  values  of  rc  (0  <  rc  <  1  minute).  Higher  concentrations  of  TCDD  are  predicted 
for  smaller  values  of  rc  due  to  the  decrease  in  the  length  of  time  before  input  of  TCDD  from 
the  liver  compartment  appears  at  the  site  of  measurement  in  the  venous  blood.  This  would 
suggest  that  models  which  ignore  the  circulatory  delay  may  result  in  an  over-prediction  of  TCDD 
concentration.  The  predicted  concentration  of  TCDD  in  the  blood  at  the  exit  from  the  liver, 
Cs(f,  1),  is  shown  in  Figure  10. 

4.3.3  Dependence  on  the  Axial  Dispersion  Number 

In  Section  4.3,  we  explained  that  the  magnitude  of  the  axial  dispersion  number  does  not  affect  the 
arterial/venous  blood  concentration,  Ca,  on  the  initial  delay  interval  [0,  rc].  On  [rc,  2 rc],  however, 
we  begin  to  see  the  effects  of  liver  concentrations  on  the  amount  of  TCDD  in  the  arterial/venous 
blood  supply.  For  T>n  very  small  (Xbv  — >  0),  the  model  is  driven  predominantly  by  convective 
forces;  that  is,  the  predictions  are  closely  related  to  that  of  a  plug-flow  model.  As  1>n  increases 
(Pjv  — >  oo),  the  model  approaches  one  with  a  well-mixed  state,  with  uniform  concentrations 
along  the  length  of  the  cylinder,  so  that  input  to  the  cylinder  at  :t:  =  0  is  immediately  distributed 
along  its  length.  This  results  in  higher  TCDD  concentrations  along  the  cylinder  as  compared  to 
those  predicted  with  lower  values  of  Pjv-  Figures  11  and  12  depict  the  variations  in  predicted 
concentrations  for  C'b  as  a  function  of  x,  after  one  minute  and  one  hour,  respectively.  The  effects 
of  these  variations  on  the  arterial  blood  compartment  are  shown  in  Figure  13,  as  Dpf  varies  from 
5  to  1000. 

4.3.4  Remarks  on  the  Magnitude  of  the  Dispersion 
Number 

Roberts  and  Rowlands  estimated  the  dispersion  number  to  be  between  0.1  and  0.2  for  red  blood 
cells  and  other  non-eliminated  solutes  following  a  bolus  input  [12]  and  0.2-0. 5  at  steady  state  fol¬ 
lowing  a  constant  input  concentration  [19].  They  noted,  however,  that  it  is  likely  that  dispersion 
numbers  obtained  for  extracted  solutes  will  generally  be  higher  than  those  for  substances  that 
are  not  eliminated  [12].  Moreover,  binding  of  solutes  to  cellular  constituents  results  in  longer 


11 


4 


redicted  concentrations  of  TC 
the  delay  period,  rc,  for  0  <  t 


CB(t,x)  for  t  =  0.01 6667  hrs,  tq  =  0.01 6667  hrs,  Nx  =  1 6 


distance,  x  (dimensionless) 


Figure  11:  Predicted  concentrations  of  TCDD  in  the  liver  blood  compartment  as  a  function  of 
distance,  x ,  and  the  axial  dispersion  number,  Vjy,  at  time  t  =  1  minute. 


CB(t,x)  for  t  =  1  hrs,  xc  =  0.01 6667  hrs,  Nx  =  1 6 


0  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9  1 

distance,  x  (dimensionless) 

Figure  12:  Predicted  concentrations  of  TCDD  in  the  liver  blood  compartment  as  a  function  of 
distance,  x,  and  the  axial  dispersion  number,  V jv,  at  time  t  —  1  hour. 
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Ca  (t):  Dependence  on  DN  for  nx  =  16,  t,  =  0.016667  hr 


Figure  13:  Predicted  concentrations  of  TCDD  in  the  combined  arterial/venous  blood  compart¬ 
ment  as  a  function  of  the  axial  dispersion  number,  'ZI/v . 


residence  times  in  the  liver;  solutes  that  are  highly  bound  in  the  cells  are  retained  in  the  liver  for 
long  periods  of  time  while  those  that  are  confined  to  the  sinusoidal  bed  pass  through  the  liver 
quickly  [12], 

The  dispersion  number  is  dependent  on  blood  flow  rates,  protein  binding,  elimination  and 
hepatic  cell  permeability,  and  is  probably  also  species-dependent  [12],  In  our  investigations  of 
the  qualitative  behavior  of  the  TCDD  model,  we  have  used  dispersion  numbers  varying  from  0.5 
to  1000.  In  simulations  where  the  dispersion  number  was  fixed  and  other  parameters  (such  as  the 
circulatory  delay,  rc)  were  varied,  we  chose  to  use  values  of  T> jv  at.  least  an  order  of  magnitude 
higher  than  Roberts  and  Rowlands’s  estimate  due  to  the  likelihood  that  their  estimates  are  too 
low  for  our  model.  TCDD  is  highly  lipophilic  [20]  and  enters  the  cells  easily,  where  it  binds  to 
two  intracellular  proteins,  the  Ah  receptor  and  CYP1A2.  The  dispersion  model  used  in  the  data 
analysis  of  Roberts  and  Rowlands  assumed  first-order  binding  to  cellular  constituents.  The  more 
complex  nonlinear  chemical  kinetics  in  our  model  may  arguably  correspond  to  a  higher  dispersion 
number.  In  addition,  the  majority  of  the  parameters  taken  from  the  literature  (e.g.,  [2,  5,  21,  22]) 
were  estimated  by  fitting  well-stirred  models  to  experimental  data.  Therefore,  although  our 
investigations  here  were  of  a  qualitative  nature,  the  use  of  a  high  dispersion  number  may  be 
quantitatively  more  in  keeping  with  those  parameter  estimates.  The  magnitude  of  the  dispersion 
number  for  the  TCDD  model  is  a  question  that  is  most  properly  addressed  in  the  context  of  a 
parameter  identification  problem,  which  would  involve  fitting  the  model  to  experimental  data. 

5  Future  Directions 

In  this  presentation  we  studied  the  qualitative  behavior  of  the  TCDD  model  and  its  dependence 
on  two  parameters,  the  axial  dispersion  number,  T>n,  and  a  circulatory  delay,  rc.  We  found 
that  predicted  concentrations  of  TCDD  in  the  liver  increased  with  either  increasing  dispersion 
number  or  decreasing  circulatory  delay.  Furthermore,  the  model  predicts  spatial  variations  in 
cellular  concentrations  of  dioxin  at  low  dispersion  numbers,  a  feature  that  the  most  commonly 
used  model,  the  well-mixed  model,  does  not  provide. 

Further  analysis  of  the  qualitative  behavior  of  the  system,  including  the  effect  of  the  induction 
delay,  is  currently  underway.  The  model  will  be  adapted  for  spatially-varying  basal  synthesis  and 
induction  of  CYP1A2,  thereby  permitting  comparisons  to  the  multi-compartment  liver  model  of 
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Andersen  et  al.  [6].  The  development  of  these  types  of  advanced,  spatially-dependent  pharma¬ 
cokinetic  and  pharmacodynamic  models  may  play  an  important  role  in  the  development  of  risk 
assessment  protocols  for  dioxin  and  other  environmental  toxins. 
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Nomenclature 


Abbr. 

Description 

Units 

Mi{t) 

exit  boundary  condition  term  resulting  from  flux  balance 

nmol /hr 

ca 

total  arterial  and  venous  blood  TCDD  concentration 

nmol/ L 

cAh 

concentration  of  available  Ah  receptor  protein 

in  hepatocytes 

nmol/ L 

C,\h  i 

concentration  of  Ah  receptor-TCDD  complex 

in  hepatocytes 

nmol/  L 

CB 

total  concentration  of  TCDD  in  liver  blood 

nmol/  L 

Cpr 

concentration  of  available  CYP1A2  in  hepatocytes 

nmol/  L 

Cpr-T 

concentration  of  CYP1A2-TCDD  complex  in  hepatocytes 

nmol/  L 

CuB 

concentration  of  unbound  TCDD  in  liver  blood 

nmol/ L 

CuH 

concentration  of  unbound  TCDD  in  hepatocytes 

nmol/ L 

VN 

axial  dispersion  number 

fuB 

fraction  of  TCDD  unbound  in  the  blood 

fuij 

fraction  of  TCDD  unbound  in  space  of  Disse 

m 

input  concentration  of  TCDD  at  time  t 

nmol/ L/h 

IpT 

maximum  rate  of  synthesis  of  CYP1A2 

in  the  presence  of  TCDD 

nmol/ L/h 

k+ 1 

association  rate  constant  of  TCDD  and  Ah  receptor 

L/nmol/h 

k+ 2 

association  rate  constant  of  TCDD  and  CYP1A2 

L/nmol/h 

k- i 

dissociation  rate  constant  of  Ah  receptor-TCDD  complex 

/hr 

k- 2 

dissociation  rate  constant  of  CYP1A2-TCDD  complex 

/hr 

k3 

apparent  first-order  metabolic  clearance  rate  of  TCDD 

/hr 

kd(Ah) 

rate  constant  for  thermal  inactivation  of  Ah  receptor  protein 

/hr 

kd(Pr ) 

rate  constant  for  degradation  of  CYP1A2 

/hr 

ke 

rest  of  body  elimination  term 

/hr 

k§  ( Ah ) 

rate  constant  for  synthesis  of  Ah  receptor  protein 

nmol/ L/h 

ks(Pr) 

rate  constant  for  basal  synthesis  of  CYP1A2 

nmol/ L/h 

P 

permeability  coefficient  of  the  hepatocytes  to  TCDD 

L  /  hr 

Q 

volumetric  flow  rate  of  liver  blood 

L  /  hr 

Qa 

volumetric  flow  rate  of  venous  blood 

L/hr 

Tc 

circulation  delay 

hr 

Tr 

lag  time  between  TCDD  binding  to  Ah  receptor  and 

cellular  response  of  CYP1A2  induction 

hr 

Va 

combined  arterial  and  venous  blood  volumes 

L 

Vb 

liver  blood  volume 

L 

Vd 

Disse  space  volume 

L 

Vh 

hepatocyt.e  volume 

L 

X 

dimensionless  spatial  variable  ( x  =  0  corresponds  to  the  inlet,  x  = 

1  to  the  outlet) 
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